{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "from htmd import *\n",
    "htmd.config(viewer='webgl')\n",
    "os.chdir('/webdata/73hboiwia98hdj209jq0/')  # Skip this command."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# System Building Trypsin-Benzamidine"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "by [Stefan Doerr](https://github.com/stefdoerr)\n",
    "&nbsp;\n",
    "\n",
    "&nbsp;\n",
    "\n",
    "**Download all the required files for the tutorial from this** [link](http://pub.htmd.org/73hboiwia98hdj209jq0/bentryp.tar.gz)\n",
    "\n",
    "You can watch the presentation here:\n",
    "\n",
    "[![](http://pub.htmd.org/73hboiwia98hdj209jq0/bentryp_youtube.png)](https://youtu.be/DF9cHKBX19A?t=5m42s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Obtain structures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "glob('bentryp/*')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "prot = Molecule('3PTB')\n",
    "prot = Molecule('bentryp/trypsin.pdb')\n",
    "prot.view()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Clean structures"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Water molecules and one calcium ion present in the PDB crystal structure were also read by HTMD. Start by removing the ligand:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "prot.filter('chain A and (protein or water or resname CA)')\n",
    "prot.view()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "prot.remove('resname BEN')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Define segments"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "To build a system in HTMD, we need to separate the chemical molecules into separate segments. This prevents the builder from accidentally bonding different chemical molecules and allows us to add caps to them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "prot.set('segid', 'P', sel='protein')\n",
    "prot.set('segid', 'W', sel='water')\n",
    "prot.set('segid', 'CA', sel='resname CA')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Center the protein to the origin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "prot.center()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Let's add a ligand!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ligand = Molecule('./bentryp/benzamidine.pdb')\n",
    "ligand.center()\n",
    "ligand.view(sel='resname MOL', style='Licorice')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But the ligand is located inside the protein...\n",
    "We would like the ligand to be:\n",
    "\n",
    "* At a certain distance from the protein\n",
    "* Rotated randomly, to provide different starting conditions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Let's reposition the ligand then"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from moleculekit.util import uniformRandomRotation\n",
    "ligand.rotateBy(uniformRandomRotation())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This took care of the ligand rotation. But now we still need to position it far from the protein.\n",
    "We need to find out the radius of the protein:\n",
    "\n",
    "![](http://pub.htmd.org/73hboiwia98hdj209jq0/maxdist.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from moleculekit.util import maxDistance\n",
    "D = maxDistance(prot, 'all')\n",
    "print(D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "D += 10\n",
    "ligand.moveBy([D, 0, 0])  # Move the ligand 10 Angstrom away from the furthest protein atom in X dimension\n",
    "ligand.rotateBy(uniformRandomRotation())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Don't forget to add a segment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ligand.set('segid','L')\n",
    "ligand.set('resname','MOL')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Mix it all together"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "mol = Molecule(name='combo')\n",
    "mol.append(prot)\n",
    "mol.append(ligand)\n",
    "mol.reps.add(sel='protein', style='NewCartoon', color='Secondary Structure')\n",
    "mol.reps.add(sel='resname MOL', style='Licorice')\n",
    "mol.view()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Solvate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "> Water is the driving force of all nature.\n",
    ">\n",
    "> -Leonardo da Vinci\n",
    "\n",
    "![](http://pub.htmd.org/73hboiwia98hdj209jq0/waterbox.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "D = D + 5\n",
    "smol = solvate(mol, minmax=[[-D, -D, -D], [D, D, D]])\n",
    "smol.reps.add(sel='water', style='Lines')\n",
    "smol.view()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Build"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "charmm.listFiles()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Build and ionize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "topos  = ['top/top_all22star_prot.rtf', './bentryp/benzamidine.rtf']\n",
    "params = ['par/par_all22star_prot.prm', './bentryp/benzamidine.prm']\n",
    "\n",
    "molbuilt = charmm.build(smol, topo=topos, param=params, outdir='/tmp/build')\n",
    "\n",
    "molbuilt.view(sel='protein', style='NewCartoon', \n",
    "              color='Secondary Structure', hold=True)\n",
    "molbuilt.view(sel='resname MOL', style='Licorice', hold=True)\n",
    "molbuilt.view(sel='ions', style='VDW', hold=True)\n",
    "molbuilt.view(sel='water', style='Lines')"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
